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The North Greenland Ice Core Project (NorthGRIP) provides paleoclimatic in- 
formation back to at about 120 kyr before present (Dahl- Jensen and others, 2002). 
Each year, precipitation on the ice sheet covers it with a new layer of snow, which 
gradually transforms into ice crystals as the layer sinks into the ice sheet. The size 
distribution of ice crystals has been measured at selected depths in the upper 880 m 
of the NorthGRIP ice core (Svensson and others, 2003b), which covers a time span of 
5300 years. The distributions change with time toward a universal curve, indicating 
a common underlying physical process in the formation of crystals. We identify this 
process as an interplay between fragmentation of the crystals and diffusion of their 
grain boundaries. The process is described by a two-parameter differential equation 
to which we obtain the exact solution. The solution is in excellent agreement with 
the experimentally observed distributions. 
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The NorthGRIP drilling location (75.10 N, 42.32 W) is situated on an ice ridge with an 
ice thickness of 3085 m. The mean annual temperature is -32 deg. Celsius and the annual 
accumulation in the area is on the average 19.5 cm ice equivalent. In the upper 80 m of 
the ice sheet, the firn, the snow gradually compacts to a close packing of ice crystals of 
typical sizes 1 to 5 mm. We apply crystal size distributions obtained from fifteen vertical 
thin sections of ice evenly distributed in the depth interval 115 - 880 m (Svensson and 
others, 2003b). The thin sections have dimensions of 20 cm x 10 cm (height x width) and 
a thickness of 0.4±0.1 mm. 




FIG. 1: An image of a 10 x 10 cm 2 vertical thin section of ice from the depth 115 m. The 
section is viewed between two crossed linear polarizers and the different colors represent almost 
3000 individual crystals with various orientations of the crystal optical c-axes. 

Digital images of ice thin sections placed between crossed linear polarizers have been 
used to map the dimensions of individual ice crystals in the sample (Fig.[IJ. The ages of 
the considered samples are all less than 5300 years (Johnsen and others, 2001), and the 
temperature of the ice in this period can be assumed constant (Dahl- Jensen and others, 
1998). 

Fig. 2 shows size distributions of ice crystals at selected depths down to 880 m (5300 years 
before present (B.P.)). The crystal size is defined as the vertical extent of a crystal, which 
is estimated as the height of the minimal and vertical aligned rectangle, which encloses the 
individual crystal. Using the horizontally measured sizes of the ice crystals from the thin 
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FIG. 2: Distributions of ice crystal sizes at depths 115m, 165m, 220m, 330m, 440m and 605m. 
The crystal size is defined as the maximum vertical extent of the individual crystals. The black 
lines are the measured histograms and the smooth lines are the temporal evolution predicted by 
eq. © starting from the initial distribution at 115m. The total counts of ice crystals decreases 
with depth (due to the overall increase of sizes) until the steady state is reached. 

sections, one obtains similar distributions. The thin sections used in determining the ice 
size distributions have typical thicknesses of 0.4±0.1 mm, and ice crystals with all length 
scales smaller than this thickness will consequently be either smoothed away when cutting 
the ice or be invisible when viewed in the linear polarizers. These small crystals will be seen 
as part of the larger crystals and not as individual crystals. This corresponds to an average 
increase in the sizes of the large crystals, i.e. the lengths we measure x are equal to the 
real length x plus some noise e of positive mean, x = x + e. We have assumed that this 
noise is sufficiently peaked around half the thickness of the thin section, such that we simply 
put e = 0.2mm. This value of e has been used in the comparisons between theoretical and 
experimental values. 

Each distribution exhibits a pronounced peak, indicating a typical crystal size at each 
depth, followed by an exponentially decaying tail of relatively large crystals. The mean size 
becomes larger with depth and thus time until it reaches a limiting crystal size (Thorsteinsson 
and others, 1997; Li and others, 1998). The most important physical process behind such 
growth is diffusion of grain boundaries between the ice crystals. Inhomogeneities on the 
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crystal boundaries, characterized by small radii of curvature, tend to be smoothened out as 
time progresses. This causes the inhomogeneities to be incorporated into the larger crystals 
leading to an overall growth of the mean crystal size (Alley and others, 1986; Paterson, 
1994). This approximative description leads to the so-called Normal Grain Growth law that 
predicts the following temporal dependence of the mean crystal size x, 

(x 2 )(t)=x 2 + kt, (1) 

where x is the initial mean crystal size and k is the crystal growth rate. This approach was 
successfully applied for the GRIP ice core by Thorsteinsson and others, 1997, who fitted 
their experimental observations well back to 2500 B.P. By definition, however, the diffusion 
does not cause the crystal size to reach a limiting size, which is in disagreement with the 
actual observations that clearly indicates a saturation in (x 2 ) after 2500 B.P. (Alley and 
Woods, 1996; Gow and others, 1997; Thorsteinsson and others, 1997). 

The important missing part of the physical mechanism in this model is the fact that 
all crystals, and in particular the large ones, are subjected to fragmentation processes or - 
as often referred to in ice physics - polygonization or rotation recrystallization (Alley and 
others, 1995; Thorsteinsson et al., 1997). We here include this process in the description 
which then leads to a physical model with a balance between the diffusion of grain boundaries 
and the fragmentation of crystals causing the mean crystal size to reach a steady state after 
2500 B.P. 

Another process that has been proposed as important for the production of small crystals 
is nucleation, or the formation of new grains, which may have c-axis orientations that differ 
from the dominating vertical orientation of the surrounding crystals. Measurements of the 
NorthGRIP fabric (c-axis orientations) suggest, however, that nucleation does not seem 
to play an important role in the Holocene (Wang and others, 2002; Svensson and others, 
2003b), and nucleation is not included in the model. 

Our proposed model is formulated as a rate equation in the quantity N(x,t), which is 
the density of ice crystals of length x at time t measured B.P. In this study x is chosen as 
the vertical extent of the crystals. At a given time, N(x, t) can be increased or decreased 
by diffusion with a diffusion constant D. It can receive fragments of size x from fragmenta- 
tions of larger crystals and it can decrease by its own fragmentation. The fragmentation is 
defined as a rate / in length and time, i.e. for a given time step dt the average number of 
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fragmentation events over a length L is fLdt. 

The fragmentation rate /, and the diffusion constant D, will depend on temperature, 
non-hydrostatic stresses in the ice and moreover be sensitive to the impurity content of the 
ice (Alley and others, 1986). 

These diffusion and fragmentation processes lead to an integral-differential equation on 
the form 

' POO 

-fxN(x,t)+2f N(x',t)dx' (2) 

J X 



dN(x,t) _ D d 2 N(x,t) 



dt dx 2 

Here, the first term on the right hand side is a diffusion term that corresponds to the 
grain boundary diffusion of the Normal Grain Growth law, i.e. if we only include this term 
in the model we reproduce for large times the parabolic behavior in (0). Note that the 
equation describes the dynamics of the full assembly of crystals and that the mean square 
of the crystal sizes therefore follows from 

poo / / poo \ 

(x 2 )(t) = j^ x 2 N(x,t)dx / (J N(x,t)dx) . (3) 

If we only include the diffusion term in eq. (J2J), we get for large times that (x 2 )(t) ~ ADt. 
The behavior at large times differs from normal diffusion by a factor of two because we use 
the absorbing boundary condition that for all times N(0,t) = 0. 

The last term in (J2J) is the contribution from fragmentation of larger crystals into crystals 
of size x. Combining the two ways a crystal of size x' > x can produce a fragment of size 
x with the assumption that there is a uniform probability, 1/x' for where the crystal break, 
we get 

fx'N{x') ■ - r (4) 

x 

where fx'N(x') is the number of crystals of size x' that fragments per time and 2/x' is 
the probability for generating a fragment of size x. Note that the assumption of a uniform 
probability already was made in the definition of the rate /. If we integrate over all 
crystals larger than x we achieve the last term in (0). 

Dividing through by / in eq. (J2J), we obtain an expression on the right hand side which 
depends on one parameter a = D / f only. The integral-differential equation has analytic so- 
lutions, .B(^7§), which are explicitly given in terms of the Bessel function Ki and eigenvalues 
A, 



x + A\ x z l 2 
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The function Ki(x) can be written as a sum of the more frequently used Bessel functions 
I(x), Ki(x) = ^=(I_i(x) — Ii(x)). The boundary condition, N(0,t) = for all times t, 
implies that only a discrete set of non-positive values for A is allowed. They are found 
by solving the equation B (^73) = 0. The largest eigenvalues are Ao/a 1 ^ 3 = 0, Ai/a 1//3 = 
-2.338 . . . , A 2 /a 1/3 = -4.088 ........ The general solution can be written as a linear com- 
bination of the eigenf unctions, 

N(x,t) = f:c n B(^-)e^ t -^ (6) 



a 



n=0 

where the c n 's can be computed from the distribution at time t = t (for further details 
see the work of (Ferkinghoff-Borg and others, 2003)). The fact that we have no positive 
eigenvalues guarantees the existence of a steady state solution, which is defined by the only 
surviving term (corresponding to Ao = 0) for high ages, 

N(x,t) ~ B(x/a 1/3 ) for t -> 00. (7) 

The characteristic time, r, of the exponential growth towards steady state is given by the 
second largest eigenvalue, Ai, i.e. 

r=-l/(Ai/)« (2.338- /-a 1 / 3 )" 1 - 

When the dynamics has reached steady state, and the mean crystal size has saturated, the 
distribution is described by the single parameter, a 1//3 , which defines the characteristic length 
scale of the system. In particular, one finds that the mean vertical size in the steady state 
is 

<1>oo = 3*/=»V3. 

x 1 r(3/2) 

In Fig. El we show the mean vertical size of the ice crystals, (x)(t) as function of time. The 
dots are the experimental values and the solid line is an exponential fit corresponding to the 
two leading terms in the solution, 

(x) 



e -(t-t )/r 



(8) 



where (x)q is the observed average length at time to = 500 years. From the figure, we 
estimate the characteristic time, r = 600±70 years and the average length in the steady state 
(^)oo = 2.9±0.1 mm. The two parameters correspond to an effective fragmentation rate and 
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a diffusion constant of respectively / = (5.2±0.6)-10~ 4 mm^'-yr 1 and D = (1.4±0.2)-10 -3 
mm yr -1 . Following eq. (3) we obtain the Normal Grain Growth crystal growth rate 
k = 4D = (5.5 ± 0.8) ■ 10~ 3 mm 2 yr _1 , which corresponds very well with the known values 
for GRIP, k = 5.6 • 10~ 3 mm 2 yr _1 (Thorsteinsson and others, 1997), and for NorthGRIP, 
k = 5.8-10~ 3 mm 2 yr _1 (Svensson and others, 2003b). In this work we are using the maximum 
vertical extent of the crystals, so that no correction value for the sectioning effect is needed. 

o 
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FIG. 3: The mean vertical size of the ice crystals shown versus their age in years B.P. The smooth 
line shows the best fit predicted from our dynamical description of ice crystal growth. From the 
fit we read off the diffusion constant, D ~ 1.4 • lCP 3 mm 2 -yr -1 , and fragmentation rate, / ~ 
5.2 • 10~ 4 yr^ 1 -mm _1 . The time scale is taken from ref. (Johnsen and others, 2001). 

Using these estimates, we can predict the time evolution of crystal sizes from any initial 
distribution. The solid lines in Fig. |21 show the time evolution of the distribution observed 
at time t = 500 years (115 m depth), in good agreement with the experimental results. 

When comparing the model to the vertical crystal sizes we should in principle take into 
account the vertical compression of crystals due to the ice flow. Had we chosen the horizontal 
crystal sizes for comparison, we should have taken into account the corresponding elongation. 
However, observations show that the crystal flattening, which is defined as crystal width 
divided by crystal height, is less than 15% for the applied samples (Svensson and others, 
2003b), so in order to keep things reasonable simple, we simply ignore this effect. 

During an annual cycle, the impurity content and the average crystal size show important 
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variations with the highest dust load and the smallest crystals appearing during spring 
(Svensson and others, 2003a). Because the size of the applied samples (20 cm depth) is close 
to the annual accumulation at NorthGRIP (19.5 cm ice equivalent), the seasonal variability 
will to a first order be averaged out in the experimental data. Still, some inter-annual 
variations in the impurity content of the ice are observed (Svensson and others, 2003a). 
These fluctuations may explain the small discrepancies between model and observations in 
Figs. 121 and El since the values of D and / are determined by the average concentration. We 
can take these fluctuations into account for the older samples, which have almost reached 
a steady state, by rescaling the crystal size distributions. In Fig. H] we demonstrate that 
this rescaling results in a universal curve for all the size distributions, described by the 
steady state Bessel function solution, which is significantly different from the widely used 
log-normal distribution. The data-collapse in the figure, i.e. the fact that all the rescaled 
distributions have the same shape, clearly supports the one-parameter nature of the steady 
state solution (jJJ). 



o 
o 
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FIG. 4: The figure shows a "data collapse" of the size distributions as a consequence of a rescaling, 
x = (log a; — (logx))/cr(log(x)), i.e. the shown distributions have zero mean and unit standard 
deviation. The lines correspond to the eight data points in Fig. [31 of the oldest samples (t > 2500 
years) and the black line on top is the steady-state solution of eq.(l). We use the rescaling x in 
order to improve the resolution around the smallest crystal sizes and note that the steady-state 
solution is transformed accordingly. 
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By formulating eq. (J2J) we thus provide a comprehensive description of the dynamical 
processes of ice crystals. Our model improves the Normal Grain Growth law based on 
diffusive growth alone and explains the reach of a steady state for the mean crystal size 
by means of fundamental physical processes. The very good agreement between data and 
model suggest that polygonization or fragmentation is the dominating process in production 
of small crystals rather than nucleation of new crystals, which is a process not included in 
the model. The suggested interplay between the diffusion and the fragmentation in the 
crystal dynamics is believed to be a central ingredient in many other systems in nature (like 
in geological and perhaps biological processes, see (Ferkinghoff-Borg and others, 2003)) and 
could provide a useful tool to predict the distribution curves of fragmented pieces. 



[1] Alley, R. B., J. H. Perepezko and C. R. Bentley. 1986. Grain growth in polar ice: II. applica- 
tion./. Glaciol, 32(112), 425-433. 

[2] Alley, R. B., A. J. Grow and D. A. Meese. 1995. Mapping c-axis fabrics to study physical 
processes in ice. J. Glaciol, ^i(137), 197-203. 

[3] Alley, R. B. and G. A. Woods. 1996. Impurity influence on normal grain growth in the GISP2 
ice core. J. Glaciol, 42(U1), 255-260. 

[4] Dahl-Jensen, D. and 6 others. 1998. Past temperatures directly from the Greenland ice sheet. 
Science, 282, 268-271. 

[5] Dahl-Jensen, D. and 8 others. 2002. The NorthGRIP deep drilling program. Ann. Glaciol, 
35, 1-4. 

[6] Li J., T. H. Jacka and V. Morgan. 1998. Crystal-size and microparticle record in the ice core 
from Dome Summit South, Law Dome, East Antarctica. Ann. Glaciol, 27, 343-348. 

[7] Ferkinghoff-Borg, J., M.H. Jensen, J. Mathiesen, P. Olesen and K. Sneppen. 2003. Competition 
between Diffusion and Fragmentation: An Important Evolutionary Process of Nature. Phys. 
Rev. Lett, 91, 266103. 

[8] Gow, A. J. and 6 others. 1997. Physical and structural properties of the Greenland Ice Sheet 
Project 2 ice core: A review. J. Geophys. Res., 102(C12), 26,559-26,575. 

[9] Johnsen, S. J. and 8 others. 2001. Oxygen isotope and palaeotemperature records from six 
Greenland ice-core stations: Camp Century, Dye-3, GRIP, GISP2, Renland and NorthGRIP. 



10 



Journal of Quaternary Science, 16(A), 299-307. 
[10] Paterson, W. S. B. 1994. The Physics of Glaciers, New York, Pergamon. 
[11] Svensson, A., P. Baadsager, A. Persson, C. S. Hvidberg and M.-L. Siggaard-Andersen. 2003a. 

Seasonal variability in ice crystal properties at NorthGRIP: A case study around 301m depth. 

Ann. Glaciol, 37. 

[12] Svensson, A. and 6 others. 2003b. Properties of ice crystals in NorthGRIP late to middle 

Holocene ice. Ann. Glaciol, 37. 
[13] Thorsteinsson, T., J. Kipfstuhl and H. Miller. 1997. Textures and fabrics in the GRIP ice core. 

J. Geophys. Res., 102(C12), 26,583-26,599. 
[14] Wang, Y., T. Thorsteinsson, J. Kipfstuhl, H. Miller, D. Dahl-Jensen and H. Shoji. 2002. A 

vertical girdle fabric in the NGRIP deep ice core, North Greenland. Ann. Glaciol, 35, 515-520. 



